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Abstract 

O^ ' We present a spectrally accurate numerical method for finding non-trivial time- 

periodic solutions of non-linear partial differential equations. The method is based 
on minimizing a functional (of the initial condition and the period) that is positive 
unless the solution is periodic, in which case it is zero. We solve an adjoint PDE to 
^ ' compute the gradient of this functional with respect to the initial condition. We include 

additional terms in the functional to specify the free parameters, which, in the case of 
the Benjamin-Ono equation, are the mean, a spatial phase, a temporal phase and the 
real part of one of the Fourier modes at i = 0. 
CN ■ We use our method to study global paths of non-trivial time-periodic solutions con- 

necting stationary and traveling waves of the Benjamin-Ono equation. As a starting 
guess for each path, we compute periodic solutions of the linearized problem by solving 
^Q , an infinite dimensional eigenvalue problem in closed form. We then use our numerical 

fT^ ' method to continue these solutions beyond the realm of linear theory until another trav- 

eling wave is reached. By experimentation with data fitting, we identify the analytical 
f~^ ' form of the solutions on the path connecting the one-hump stationary solution to the 

00 ! two-hump traveling wave. We then derive exact formulas for these solutions by explic- 

^D ■ itly solving the system of DDEs governing the evolution of solitons using the ansatz 

suggested by the numerical simulations. 
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1 Introduction 

A fundamental problem in the theory of ordinary and partial differential equations is to 
determine whether the equation possesses time-periodic solutions. Famous examples of 
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ordinary differential equations witli periodic solutions include the Brusselator |FB85l[Gov00[ 
IStrOO) and the three-body problem |Are63[ IKNWOO] . In partial differential equations, time- 
periodic solutions can be "trivial" stationary or traveling waves, or can be genuinely time- 
periodic. Such a problem can be studied in either the forced or unforced context. Forced 
problems include an external force in the PDE that is time-periodic; solutions with the 
same period are then sought. In the unforced problem, the period is one of the unknowns. 
In this work, we present a numerical method for finding genuinely time-periodic solutions of 
the unforced Benjamin-Ono equation with periodic boundary conditions. These solutions 
have many remarkable properties, which we will describe. 

Our work is motivated by the calculations of Hou, Lowengrub, and Shelley for the vortex 
sheet with surface tension |HLS941 IHLS97J . and by the analysis of Plotnikov, Toland and 
looss |PT0HITPT05| for the water wave. Hou, Lowengrub, and Shelley developed an efficient 
numerical method to solve the initial value problem for the vortex sheet with surface tension. 
They performed calculations for a variety of initial conditions and values of the surface 
tension parameter, and found many situations in which the solutions appear to be close to 
time-periodic. They did not, however, try to measure the deviation from time-periodicity 
or attempt to vary the initial conditions to reduce this deviation. Plotnikov, Toland, and 
looss have proved the existence of time-periodic water waves, without surface tension, in 
the case of either finite or infinite depth. This is proved using a version of the Nash-Moser 
implicit function theorem. Their work includes no computation of the water waves. We aim 
to get a firmer handle on these solutions with an explicit calculation. To this end, in the 
present work, we develop a general numerical method for finding time-periodic solutions of 
nonlinear systems of partial differential equations and eventually plan to use this method 
for the vortex sheet and water wave problems. In fact, during the review process of the 
present work, we have succeeded in computing several families of time-periodic solutions of 
the vortex sheet with surface tension |AW09] . 

The Benjamin-Ono equation, developed in |Ben671 [DA67t [Ono75) . is a model equation 
for the evolution of waves on deep water. It is a widely-studied dispersive equation, and 
much is known about solutions. It would be impossible to mention all results on Benjamin- 
Ono, but we mention, for example, that weak solutions exist for uq € L^ |Sau791 IGV91J . 
and that the solution exists for all time if uq E H^ |Tao04) . We chose the Benjamin-Ono 
equation as a first application for our numerical method because it is much less expensive 
to evolve than the vortex sheet or water wave, yet has many features in common with them, 
such as non-locality (due to the Hilbert transform in the former case and the Birkhoff-Rott 
integral in the latter two cases.) In our numerical simulations, we find that the Benjamin- 
Ono equation has a rich family of non-trivial time-periodic solutions that act as rungs in 
a ladder connecting traveling waves with different speeds and wavelengths by creating or 
annihilating oscillatory humps that grow or shrink in amplitude until they become part of 
the stationary or traveling wave on the other side of the rung. The dynamics of these non- 
trivial solutions are often very interesting, sometimes resembling a low amplitude traveling 
wave superimposed on a larger carrier signal, and other times looking like a collection of 
interacting solitons that pass through each other or bounce off each other, depending on 
their relative amplitudes. 

By fitting our numerical data, we have determined that all the solutions we have com- 
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puted are special cases of the multi-phase solutions studied by Case |Cas79| . Satsuma and 
Ishiniori |SI79| . Dobrokhotov and Krichever [DK91| . and Matsuno |Mat04j . but with special 
initial conditions (that yield periodic orbits) and a modified mean (to change their speeds 
and allow bifurcations to occur between different levels of the hierarchy of multi-phase so- 
lutions). We did not take advantage of this structure when we developed our numerical 
method; hence, our approach can also be used for non-integrable problems. We also note 
that bifurcation within the family of multi-phase solutions has not previously been dis- 
cussed in the literature, nor has the remarkable dynamics of the Fourier coefficients of these 
solutions beyond the original derivation by Benjamin [Ben67j of the form of the traveling 
waves for this equation. 

We are aware of very few works on the existence of time-periodic solutions for water 
wave model equations. Crannell has demonstrated |Cra96j the existence of periodic, non- 
traveling, weak solutions of the Boussinesq equation using a generalization of the mountain 
pass lemma of Rabinowitz. Chen and looss have proved existence of time-periodic solutions 
in a two-way Boussinesq-type water wave model |CI05] . As in |PT01] and [IPT05| . there is 
no computation of the solution in either of these studies. Cabral and Rosa have recently 
discovered a period-doubling cascade of periodic solutions for a damped and forced version 
of the Korteweg-de Vries equation |CR04] . They use a Fourier pseudospectral method for 
the spatial discretization and a first order semi-implicit scheme in time. To find periodic 
solutions, they use a secant method on a numerical Poincare map. Whereas our approach 
is based on minimizing a functional that measures deviation from periodicity, they rely on 
the stability of the orbit to converge to a periodic solution. They stop when they find 
a solution that returns to within one percent of its initial state, whereas we resolve our 
periodic solutions to 13-15 digits of accuracy, which allows us to study the analytic form of 
the solutions. 

Water waves aside, many authors have investigated time-periodic solutions of other par- 
tial differential equations both numerically and analytically. For instance. Smiley proves 
existence of time-periodic solutions of a nonlinear wave equation on an unbounded domain 
|Smi89| ; he also develops a numerical method for the same problem [Smi90j . On a finite do- 
main, Brezis uses duality principles to prove the existence of periodic solutions of nonlinear 
vibrating strings in both the forced and unforced setting; see |Bre83| . Mawhin has written 
a survey article on periodic solutions of semilinear wave equations [Maw95j , which includes 
many references. Pao has developed a numerical method for the solution of time-periodic 
parabolic boundary- value problems |Pao01j . Pao gives various iterative schemes, but unlike 
the present work, these are not based on variational principles or the dual system. Brown 
et. al. |BKOR92] have used the popular software package AUTO tD KK91a] to study Hopf 
bifurcation and loss of stability for traveling waves of the regularized Kuramoto-Sivashinsky 
equation. They observe modulated traveling waves similar to many of the solutions we have 
found for Benjamin-Ono. And of course, time-periodic solutions of systems of ordinary 
differential equations have also been widely studied, e.g. in |Rab78[ [Rab82[ [Zeh83l IDui84] . 

The most popular methods for the numerical solution of boundary value problems gov- 
erned by ODEs, which include finding time-periodic solutions as a special case, are orthog- 
onal collocation methods |DKK91b] and shooting/multi-shooting methods |SB02j . These 
methods can also be used for PDEs via the method of lines, i.e. by discretizing space first to 
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obtain an ODE. Orthogonal collocation methods (such as implemented in AUTO) require 
a nonlinear system of equations to be solved involving all the degrees of freedom at every 
timestep. In many of our simulations, the solution at every timestep cannot even be stored 
simultaneously in memory; hence, it would be impossible to solve a system of equations for 
all these variables. For PDEs, shooting methods are also very expensive as the Jacobian of 
the functional measuring errors in the boundary conditions must be computed with respect 
to variation of the N "unknown" initial conditions, which is N times more expensive than 
computing the functional itself. 

Lust and Roose |LR96j propose an interesting solution to this problem in which a shoot- 
ing method is used for some of the degrees of freedom while iteration on a Poincare map 
is used for the other degrees of freedom. They aim to exploit the fact that many high 
dimensional systems actually exhibit low-dimensional dynamical behavior. In the present 
work, we offer an alternative strategy based on the following idea: while the Jacobian of the 
residual measuring error in each of the boundary conditions is expensive to compute, the 
gradient of the sum of squares of this residual can be computed by solving a single adjoint 
PDE. Although the full Jacobian gives more information, it is N times more expensive to 
compute. We find that it is more efficient to use a quasi-Newton method with only the 
gradient information than it is to use a full Newton method with the Jacobian. This ap- 
proach is not restricted to finding time-periodic solutions; it will work for any boundary 
value problem. Also, although we have not tried it, our method can be easily adapted to 
incorporate the idea behind multi-shooting methods, namely that by introducing "interior" 
initial conditions and boundary conditions, we can avoid a great deal of the ill conditioning 
due to nearby trajectories diverging exponentially in time. 

The closest numerical method to our own that we have found is due to Bristeau, Glowin- 
ski and Periaux |BGP98] . who developed a least squares shooting method for numerical 
computation of time-periodic solutions of linear dynamical systems with applications in 
scattering phenomena in two and three dimensions; see also |GR06| . These authors employ 
methods of control theory to compute variational derivatives, and although they only apply 
their methods to linear problems, they mention that their techniques will also work on non- 
linear problems. Our method can be considered an extension of their approach that focuses 
on the difficulties that arise due to non-linearity. In particular, we replace their conjugate 
gradient solver with a black-box minimization algorithm, (the BEGS method |NW99j ). and 
include an additional penalty function to prescribe the values of the free parameters that 
describe the manifold of non-trivial time-periodic solutions. Without this penalty function, 
the basic method is only found to produce constant solutions and traveling waves. 

This paper is organized as follows: In Section [21 we discuss spatially periodic stationary 
and traveling solutions of the Benjamin-Ono equation, the bifurcations from constant so- 
lutions to traveling waves, and the pole dynamics of meromorphic solutions. In Section [3l 
we investigate time-periodic solutions of the linearized Benjamin-Ono equation; this is the 
linearization about the stationary solutions discussed previously. To analyze the linearized 
problem, we compute (numerically) the spectrum and eigenfunctions of the relevant linear 
operator and deduce their analytic form by trial and error; the resulting formulas can be 
verified rigorously (but we omit details). 

In Section HI we describe our numerical method, which involves minimizing a non- 
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negative functional that is zero if and only if the solution is periodic. We solve an adjoint 
PDE to compute the variational derivative of this functional with respect to perturbation of 
the initial condition and use the BFGS minimization algorithm to minimize the functional. 
The Benjamin-Ono and adjoint equations are solved with a pseudo-spectral collocation 
method using a fourth order, semi-implicit Runge-Kutta scheme. We use a penalty function 
to rule out constant solutions and traveling waves, and to prescribe the free parameters of 
the manifold of non-trivial solutions. We then vary a bifurcation parameter to study the 
global properties of these non-trivial solutions. In the present work, we apply this method 
only to the Benjamin-Ono equation, but we are confident that this method is applicable to 
virtually any system of partial differential equations that possesses time-periodic solutions. 
In Section [5l we use our method to study the global behavior of non-trivial time-periodic 
solutions far beyond the realm of validity of the linearization about stationary and traveling 
waves. We will follow one such path to discover that the one- hump stationary solution is 
connected to the two- hump traveling wave by a path of non-trivial time periodic solutions. 
In Section [U we re- formulate the ODE governing the evolution of poles to reveal an exact 
formula for the solutions on the path studied numerically in Section [5j Thus, unexpectedly, 
we have proved that non-trivial time-periodic solutions bifurcate from stationary solutions 
by exhibiting a family of them explicitly. In a follow-up paper [AWj . we will classify all 
bifurcations from traveling waves, study the paths of non-trivial solutions connecting several 
of them, propose a conjecture explaining how they all fit together, and describe their analytic 
form to the extent that we are able. We end with a few concluding remarks in Section [71 

2 Stationary, Traveling and Soliton Solutions 

We consider the spatially periodic Benjamin-Ono equation with the following sign conven- 
tion: 

Ut = Huxx-uu^. (1) 

Of course, the operator H is the Hilbert transform. Recall that the symbol of H is H{k) = 
—isgD.{k). Many authors include a factor of two in the convection term and place a minus 
sign in front of H. The former change causes solutions to be scaled by a factor of 1/2 while 
the latter change has no effect as H is defined with the opposite sign with this convention. 
This equation possesses a two-parameter family of stationary solutions, namely 



1-3/3^ 4/3 [cos (x -e)-P] 

1-/32 ^ l + /32-2/3cos(x-0) ' 



u{x)='-^^+ ,:zr::i'~/\. ^ (-i</3<i,^gm). (2) 



These solutions have mean a, related to /3 via 

aW) = y^^, l/3P = j^. (3) 

Changing the sign of /3 is equivalent to the phase shift 9 ^ 9 — n. It is convenient to 
complexify (3 and define U/j to be the mean-zero part of ([2]) with /3 -^ |/3j, 9 — )• arg/3: 

stationary solution = a(/3) + u^(x), /3 = |/3|e"*^ e A = {z e C : \z\ < 1}. (4) 
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Stationary Solutions 



a = -0.8225 



a =1 -(3k/40r, k=0..18 



a = 0.544375 




27C 



Figure 1: Stationary solutions of the Benjaniin-Ono equation. 



Note that the subscript /3 does not indicate a derivative here. Several stationary solutions 
with /3 real and negative are shown in Figure [TJ The Fourier representation of n^ is simply 



■W/3,fc 




(5) 



where /3 is the complex conjugate of /3. These functions u^(x) are the building blocks for 
the meromorphic solutions discussed below. 

Note that the constant solution u = oq is also a stationary solution, as are the re-scaled 
solutions 

unAx) = Na{p) + Nup{Nx), (/3 G A, iV = 1, 2, 3, . . . ), 

which have mean oq = Na{f3). If we restrict attention to stationary solutions with even 
symmetry (i.e. with f3 real), we find that there is a pitchfork bifurcation at each positive 
integer (using the mean as a bifurcation parameter). As oq changes from N^ to N~, the 
constant solution splits, yielding two additional (A^-hump stationary) solutions, namely 
UN,p{x) with /3 = 0^. The pitchfork would be obtained by plotting the real part of the 
Ath Fourier mode versus the mean, where we observe that the Fourier representation of 
u = njv,^ (for any /3 G A) is given by 



Uk 



Aa(/3), 


A; = 0, 


2NP^I^, 


k G AZ, A; > 0, 


2A^I^I/^, 


k G AZ, A; < 0, 





otherwise. 



(6) 
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If we do not restrict attention to even solutions, the phase shift 6 acts as a second parameter 
connecting the two outer branches of the pitchfork into a two-dimensional, bowl-shaped 
sheet (plotting the real and imaginary parts of the A'^th Fourier mode versus the mean) . 

In the bifurcation problem just described, we varied the mean ao and found bifurcations 
from constant solutions to stationary solutions at the positive integers. The remainder of 
this paper deals with bifurcation from these stationary solutions to non-trivial time-periodic 
solutions and their global continuation beyond the realm of linear theory. Rather than 
varying the mean, we will hold ao G K constant and use another quantity (such as the 
period T or the real part of one of the Fourier modes of u at t = 0) as the bifurcation 
parameter. As a first step, let us consider bifurcation from constant solutions to traveling 
waves holding ao S M constant and varying T. 

All traveling wave solutions of the Benjamin-Ono equation can be found by applying 
a simple transformation to a stationary solution, and vice versa. Indeed, if u{x, t) is any 
solution of ([I]), then 

U{x,t) =u{x-ct,t) + c (7) 

is also a solution; thus, adding a constant c to a stationary solution causes it to travel to 
the right with speed c. We can parametrize these A^-hump traveling waves by their mean 
ao € M and decay/phase parameter /3 G A: 

Uao,N,l3ix,t) =UN,I3{X-Ct)+C, (c = Oq - iVa(/3)) . (8) 

If we express the period T = 27r/(A^|c|) in terms of /3 and solve for /?, we find that we can 
bifurcate from any constant solution n = ao to an A^-hump traveling solution with the same 
mean. If ao < A^, a pitchfork from the constant solution occurs at To = 2it/[N{N — ao)]; 
as we increase T from To to oo, a = [2tt/{NT) + ao]/N decreases from 1 to oq/N and 
1/3 1 varies from to -y/(l — ao/N)/{3 — ao/N). Similarly, if ao > A^, a pitchfork occurs at 
To = 27r/[A^(ao — A^)]; as we decrease T from To to 0, a = [ao — 27r/{NT)]/N decreases 
from 1 to — oo, and |/3| varies from to 1. And if ao = A^, the situation is qualitatively 
similar to the latter case, but the pitchfork occurs at To = oo, i.e. all three solutions (with 
/3 real) exist for any period T > 0. 

We remark that if the traveling waves described above have zero mean, they are a special 
case of the meromorphic A^-particle "periodic soliton" solutions described in [Gas 78) . namely 

u{x,t) - 




Z-^ ^i\x+t-xi{t)\ _ -y 



where lm{x;(0)} > and the xi(t) satisfy the system of differential equations 

-^=y , \ + y , , , (l</<Ar). (9) 

dt ^-^ Q-t{Xm-Xl) _ I Z-^ g-t{xi-Xm) _ l' \ — — J \ J 

m=l 771=1 

In our notation, we write 



xi{t) = Hog (3i{t) = 9i{t)-nogmt)\, (A = |A|e" 
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and generalize to the case that the mean oq can be non-zero. We find from ([9]) and ([7]) that 

N 

u{x, t) = ao + Y^ ^ft(t) (^) (10) 

1=1 

is a solution of ([T]) if the variables /3; € A satisfy 

A=E^^ + E^^ + ^(2^-i-"o)A, (l</<iV). (11) 

-I Pi Pm -, p; — Pm 



m^ 



The A^-hump traveling wave then has the representation 

N 

Uao,N,p{^,t) = aQ + ^up^it){x), j5i{t) = ^/pe~^''\ c = ao-Na{P), 
1=1 

where each (3i is assigned a distinct A^th root of /?. As we are interested in developing 
numerical methods that generalize to more complicated systems such as the vortex sheet 
with surface tension and the water wave, we do not exploit the existence of meromorphic 
solutions in our numerical method; however, the non-trivial time-periodic solutions we find 
do turn out to be of this form; see Section [6j 

3 Linear Theory 

We formulate the problem of finding time-periodic solutions of the Benjamin-Ono equation 
as that of finding an initial condition uq and period T such that F{uq,T) = 0, where 
F : H^ xR^ H^ IS given by 

F{uo,T) =u{-,T) -Uq, Ut = HUxx-UUx, u{-,0)=uo. (12) 

Clearly, stationary solutions are periodic with any period T. In this section, we linearize 
F about these stationary solutions and use solutions of the linearized problem as initial 
search directions to find time-periodic solutions of the nonlinear problem. Bifurcation from 
traveling waves can be reduced to this case by adding an appropriate constant and requiring 
that the period of the perturbation coincide with the period of the traveling wave (although 
there may be a phase shift involved as well). We present a detailed analysis of the traveling 
case in [AW] . 

3.1 Linearization About Stationary Solutions 

Let u = un^p be an arbitrary A^-hump stationary solution. If u{x) + u(x, t) is to satisfy ([T]) 
to first order in v, then v should satisfy 

vt = Hvxx - iuv)x- (13) 

(The exact solution satisfies vt = Hvxx — {uv)x — vvx-) Equation (J13p can be written 

Vt = iBAv, (14) 
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where the (unbounded, self-adjoint) operators A and B on H^ are defined as 

1 



A = Hd^ 



u, 



B 



rdrr.. 



To solve (|14p . we are interested in the eigenvalue problem 

BAz 



coz, 



(15) 



(16) 



so that if BA has a complete set of eigenvectors, the general solution of (J14p will be a 
superposition of functions of the form 



v{x,t) =Re{Cz{x)e''^'} 



C £ 



Of course, the eigenvalues of a composition of Hermitian operators need not be real, but for 
A and B in psp . we can compute all the eigenvalues explicitly, and they are indeed real. 
We do this numerically (which surprisingly leads us to formulas we can check analytically) 
by truncating the Fourier representations of A and B and computing the eigenvalues of the 
matrix BA. More precisely, we choose a cutoff frequency K (e.g. K = 240) and define the 
{2K - 1) X {2K - 1) matrices 



A 



kl 



\k\^kl - Uk-l = \k\Ski - ui-k, 



Bki = k6ki, 



-K<k,l<K), (17) 



where Uk was given in ^ and (5^^ = 1 if A; = / and otherwise. By carefully studying the 
eigenvalues for different values of N and /3 = —\/{l — a)/(3 — a) with a < 1, we determined 
that 

-UJN-n n < 0, 

n = 0, 

{n){N-n) l<n<N-l, 

{n + l-N){n + l + N{l-a)) n>N. 

With this numbering, the first A^ — 1 non-zero eigenvalues are independent of a: 



<^N,n 



(18) 



N 



n 



WW, 




(19) 



Note that ci;jv,iv = (2 — a)N + 1 > A^ + 1 and UN^n is strictly increasing in n for n > 
N, but uJiy,N could be less than i^N.lN/2] when A^ > 6 (and some of the eigenvalues can 
coalesce, increasing their multiplicity). Nevertheless, the ordering of the eigenvalues in p8p 
is more convenient than the monotonic ordering due to the fact that a pathway of non- 
trivial solutions connecting an A^-hump traveling wave to an A^'-hump traveling wave with 
N < N' seems to involve uj]sf „ and u^i „/ with n > N and n' < N' satisfying A^' = n + 1 
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and n' = N' — N (see |AWj .) These global reconnection formulas would be much more 
complicated if the eigenvalues were ordered monotonically. 

The zero eigenvalue ujn,o = has geometric multiplicity two and algebraic multiplicity 
three. The fact that the dimension of the kernel is independent of a indicates that there 
are no special values of the mean Na at which these A^-hump stationary solutions bifurcate 
to more complicated stationary solutions. The two eigenfunctions in the kernel of BA are 



&ix) = -Mx) 



d_ 
dd 



e=o 



(2) 
UN,I3{X - 0), zy^Q{x) 



d 



m 



UN,l3{x), 



(20) 



which correspond to translating the stationary solution by a phase or decreasing its mean, 
Na = A^(l — 3|/3p)/(l — |/3|^). There is also a Jordan chain |Wil07aj of length two associated 
with zJy'q (x), namely 



(1.1)/ \ _, 
The corresponding solution of (I14p is 



BAz 



(1,1) 

N,0 



,(1,1) 



(1,0), 



v{x,t) = -izljy{x) + tz}^y{x) = 1 - tn^(x) 



Ji,o) 



d_ 

de 



(21) 



\u[x 



et)+e\, 



e=0 



i.e. this linear growth mode arises due to the fact that adding a constant to a stationary 
solution causes it to travel. The multiple eigenvalues uJN,n = ^N,N-n with 1 < n < A^ — 1 
pose a minor obstacle to obtaining explicit formulas for the eigenvectors. We eventually 
realized that because the shift operator 



Sez{x) 



zix 



9), Se., 



kl 



e-'^^Hku 



2'k/N 



(22) 



commutes with BA, the eigenspaces of BA are invariant under the action of Sg. Thus we 
can impose the additional requirement that if z is an eigenvector of BA corresponding to a 
multiple eigenvalue, then z should also satisfy 



Zfc 7^ and zi ^ Q 



k-le NIs, 



(23) 



i.e. the non-zero Fourier coefficients are equally spaced with stride length A^. Using this 
condition to make the eigenvectors unique up to scaling, we were able to recognize the 
patterns that emerge in the numerical eigenvectors (with the exception of the coefficient C 
and the j = case when n > N, which we determined analytically): 



ZN,n,k 



i + mhii)^m-^ 



k=n+jN 



C 1 + 



M\ ^i+i 



j<0 
j>0 



1 <n< A^-l 

/^ —nN 



{N-n)[n+{N-n)\l3\'2\ I ' 



ZN,n,k 



k=n-N+l+jN 



-P 



(1-I/3P) 



T^- 



1 



1 



N 
n+1 



1 + 



N{3-1) \ m-X 



n.+l 



/3^ 



i 


<0 ^ 


3 


= 


j 


> 1 



{n>N). (24) 
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These formulas can be summed to obtain ZN,n{x) as a rational function of e*^, but we prefer 
to work with the Fourier coefficients. Note that as n ^^ oo (holding A^ fixed), the index 
fe = n — A^ + lof the first non-zero Fourier mode increases to infinity. The eigenvectors 
corresponding to negative eigenvalues ojn-u with n > 1 satisfy Zf^^-n{x) = z^^nix), so 
the Fourier coefficients appear in reverse order, conjugated: ZN,-n,k = Z]y^n,-k- When /3 
is real, the Fourier coefficients are real and zj\[^-nix) = ziy,n{—x). We have verified the 
formulas (llSp and (I24p analytically, and can also prove that the Fourier representation of 
these eigenvectors (together with the associated vector corresponding to the Jordan chain) 
form a Riesz basis for ^^(Z); hence, we have not missed any eigenvalues. 

3.2 Bifurcation from Stationary Solutions 

Now that we have solved the eigenvalue problem for BA, we can compute the derivative 
of the operator F in (I12p above. We continue to assume that u is an N-hump stationary 
solution so that DF = {DiF, D2F) : H^ xR^ H'^ satisfies 



DiF{u,T)vo = -^\^^^F{u + evo,T) = v{;T) - vq = [e^^^^ - /] vq, 

D2F{u,T)t = — F{u,T + eT) = 0. 
oe e=0 



(25) 



Note that vq € kei DiF{u,T) iff the solution v{x,t) of the linearized problem is periodic 
with period T. These solutions of the linearized problem serve as initial search directions in 
which to find periodic solutions of the non-linear problem. Since D2F = in the stationary 
case, the nuhspace TV of DF = {DiF, D2F) is of the form TV = TVi x M where TVi = ker DiF. 
A basis for TVi consists of the functions 

vo{x) = ReizN^nix)} and vq{x) = lm.{zN,n{x)], (26) 

where n ranges over all integers such that 

ujN,nT G 27rZ. (27) 



Negative values of n have already been accounted for in ([26]) using ZN-n{x) = ZN,n{x)- The 
n = case always yields two vectors in the kernel, namely those in (I20p . These directions 
do not cause bifurcations as they lead to other stationary solutions in the two parameter 
family. Thus, the periods at which bifurcations are expected are 

2irm 
TN,n,m = , (m,n>l). (28) 

Note that this set is dense on the positive real line since u}N,n — >• co as n — )■ 00. This 
leads to a small divisor problem when trying to apply the Liapunov-Schmidt reduction 
|GS85] IKie04j to rigorously analyze bifurcations in the set of solutions of the equation 
F{uo,T) = 0. For other problems, small divisors have been dealt with successfully using 
Nash-Moser theory; see e.g. [NirOU IPTOH llPTOSj . As most of the bifurcations responsible 
for the small divisor problem involve high frequency perturbations that are smoothed out 
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Amplitude of traveling wave versus period 




47112.6 



13.1 



Figure 2: First seven bifurcations from the constant solution u{x) = qq to traveling waves 
with n-humps. The period shown is T = 2TTm./[n{na — ao)] with m = 2n{n — 1/2). We used 
this to solve for |/3| in terms of T via ([3]). The amplitude shown is the difference between 
the maximum and minimum values of the solution, i.e. 8n|/3|/(l — |/3p). 



by the numerical discretization, we have not found small divisors to cause major difficulties 
in our ability to track paths of time-periodic solutions. 

In our numerical studies, we have found that each of the eigenfunctions ZAr_„(x) with 
n > 1 gives a direction along which a sheet of non-trivial solutions bifurcates from the N- 
hump stationary solution. More precisely, if we let a parameter e — > 0, there appear to be 
time-periodic solutions of the nonlinear problem with period T^ = TN,n,m + 0{e'^) that agree 
with the real and imaginary parts of ez^^ni^) at t = and t = —Tf,/[Am), respectively, 
to 0{e'^). This is interesting because in cases that several (even infinitely many) of the 
eigenfunctions zi^^n{x) lead to solutions of the linearized problem with the same period, 
most linear combinations of these solutions will not give a bifurcation direction. In other 
words, the eigenfunctions ZN^n{x) are already diagonalized with respect to the nonlinear 
effects of the bifurcation. For this to work out, it is essential that when ci^Af^n is a multiple 
eigenvalue, ZN,n{x) is chosen as in (f23|) to simultaneously diagonalize the shift operator 

'S'27r/Ar- 

We observe an analogous phenomenon when bifurcating from constant solutions to trav- 
eling waves; see Figure [2j When n = oq is a constant function in (|15p above, A and B are 
both diagonal matrices, the eigenvalues and eigenvectors of BA are given by 

ujn = n{\n\ - Qo), Zn{x) = e^""", (n G Z), 

and the bifurcation times are given by 



27r?n,/(n|ao — n.|), {n,m > 1). 
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Note that in this simphfied problem, the bifurcation index n turns out to be the number 
of humps. If ao = 1/2 and T = 4tt, then uj^T G 27rZ for every n, i.e. the kernel A^i 
of DiF{u,T) is the whole space H^ . Nevertheless, the traveling solutions that emerge 
from this bifurcation are no different than if ao were irrational and Mi were spanned by 
{1, cosnx,sinnx} for some fixed n — they all just happen to join together at T = 47r 
for this value of the mean. More specifically, the n-hump traveling solutions Ua^^n^p{x,t) 
defined in ([8]) above have the property that as /? — >■ (and hence a — >■ 1), a multiple m of 
their shortest period 27r/[n(na — oq)] converges to 47r. So while cosx and cos2j; are each 
bifurcation directions for initial conditions that lead to periodic solutions of the nonlinear 
problem, (cosx + cos2x) is not. 

4 The Method 

In order to compute non-trivial time periodic solutions, we define the functional 

Gtot(uo,r) = G(no,T) + ^(uo,r) (29) 

with 

G(no, T) = ]- j [u{x, T) - UQ{x)f dx (30) 

^ JO 
and look for minimizers of Gtot with the hope that the minimum value will be zero. Here 
if{uo,T) is a non-negative penalty function designed to eliminate spatial and temporal phase 
shifts and specify the mean ao and the amplitude of one of the Fourier modes at i = 0. 
Our first goal is to find an efficient method of computing the variational derivative of G. As 
usual in optimal control problems |Pir84j . there is an adjoint PDE that allows us to compute 
■^ in as little time as it takes to compute G itself. We will then use a spectral method in 
space and a fourth order semi-implicit Runge-Kutta method |CS831 IKC031 IWilOTbj in time 
to solve the Benjamin-Ono and adjoint equations to compute G, -^ and ^ in the inner 
loop of the BFGS minimization algorithm [Bro70t INW99| . 

4.1 Variational Derivative of G 

Let uo be any function in H^ (not necessarily leading to a periodic solution). Evidently, 

—G{uo,T)= [u{x,T) -uo{x)]ut{x,T)dx. (31) 

Let vq € H^ be given and define G = DiG{uq,T)vo, i.e. 



de 



d '■2- 



G{uo + £Vq,T)= / [u{x,T) - uo{x)][v{x,T) - vo{x)]dx. (32) 

£=0 Jo 



Here v{x,t) = it{x,t) = ^| u{x,t,e) with u{x,t,e) the solution of Benjamin-Ono with 
initial condition u{x,0,e) = uo{x) -\-£Vq{x). We can compute v by solving the variational 
equation 

vt = Hvxx- {uv)x, v{x,0) =voix), (33) 
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which is hnear but non-autonomous (as u depends on time in general). Our next task is to 
eUminate v{x,T) from (|32p and represent G as an inner product: 

G= - — {x)vo{x)dx. 

Jo ouo 

The idea is to define a function w{x, s) going backward in time (with s = T — t) such that 



w{x, 0) = wo{x) = u{x, T) — uq{x) 



(34) 

that 
(35) 



and then determine how w should evolve so that 



2-K /•2tt 

w{x,0)v{x,T) dx = / w{x,T)v{x,0) dx. 
Jo 



(36) 



Let us define the solution operator V{t2,ti) : H^ -^ H^ for the linearized equation (j33l) as 
the mapping that evolves an initial condition specified at time ti to the solution at time t2- 
These operators satisfy a non- autonomous, time reversible version of familiar semigroup 
properties: 

V{ti,ti)=I, V{t3,ti)=V{t3,t2)V{t2,ti), (ti,t2,t3GK). (37) 

Equation (j36p may now be written 

{wo, V{T, 0)vo) = {W{T, 0)wo, vo) (38) 

where (•, •) is the L^ inner product and we define W{s2,si) = ^(^1,^2)* with tj = T — Sj. 
It follows from ([371) that W{si,si) = I and W{s3, si) = W{s3, S2)W{s2, si). What remains 
is to determine how this non-autonomous semigroup W is generated. Taking the inner 
product of vt with w, we have 



vt{x,t)w{x,s) dx = lim 



Vit + h,t)-V{t,t) 



h 



v{x, t) ) w{x, s — h) dx 



lim / v(x,t) 
h^o ' 



W{s,s-h) 
h 



w{x,s — h)\ dx = v{x,t)ws{x,s) dx. (39) 



We learn that 

vw^ dx 



/ vtw dx = [-ffwxx - {uv)x]w dx = v[-Hwxx + uw^] dx, (40) 

i.e. w should solve the adjoint equation to (|33]) . namely 

Ws{x, s) = -Hwxxix, s) + u{x, T - s)wx{x, s). (41) 

The time reversal in the inhomogeneous term u{x, T — s) is significant. Combining this with 
([32|) and (p^ . we conclude that 



- — (x) = w{x,T) - wo{x) 
duo 



(42) 



where w solves ()4ip with initial condition (j35p . 

Remark: We emphasize that the steps we have just followed for the Benjamin-Ono equation 

can in principle be carried out for any PDE. These steps are simply: 
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1. Find the variational equation analogous to (J13p 

2. Find the appropriate adjoint equation, accounting for time-reversal. 

The details of the initial condition of the adjoint problem and the formula for ^^ depend 
on the particular functional G we choose, but they are usually straightforward to work out. 
For example, as another variant, we could define 

G{uo,T) = l [^[u{x,T/2)-u{27r-x,T/2)fdx (43) 

^ Jo 

to impose even symmetry at the half-way point. (Recall that if uq is symmetric, then 
u{2tt — x,T/2) = u{x, —T/2)). In this case we find that 

-^(x) = 2w(x,T/2), wo(x) = u(x,T/2) -u(2tt-x,T/2), (44) 

duo 

or, since vq is assumed symmetric in this formulation, j^{x) = w{x, T/2) + w{27r — x, T/2). 
In subsequent work, we will apply the methods of this paper to the vortex sheet with surface 
tension and to the water wave. Although step 2 usually amounts to a simple integration by 
parts as was done in (|40p above, the adjoint calculation can be quite involved for systems of 
PDEs with more complicated nonlinearities such as the Birkhoff-Rott integral in the vortex 
sheet problem (see |AW09] ). 

4.2 The Numerical Method 

We minimize Gtot using the BFGS algorithm |NW 99] . which is a quasi- Newton line search 
method that builds an approximate Hessian incrementally from the history of gradients it 
has evaluated. As a black box unconstrained minimization algorithm, it requires only an 
initial guess and subroutines to compute Gtot(Q) and VqGtot(9), where g G M'^ contains 
the numerical degrees of freedom used to represent uq and T. We use a solution of the 
linearized problem for the initial guess near a stationary solution or traveling wave, and 
then use linear extrapolation (or the result of the previous iteration) for the initial guess in 
subsequent calculations as we vary the bifurcation parameter. 

In our implementation, we wrote a C++ wrapper around J. Nocedal's L-BFGS Fortran 
code released in 1990, but we turn off the limited memory aspect of the code since computing 
G takes more time than the linear algebra associated with updating the full Hessian matrix. 
We do find that the algorithm converges quadratically once it gets close to a minimizer. 
Our code also makes use of the FFTW and LAPACK libraries, but was otherwise written 
from scratch. 

We represent u{x, t) spectrally as a sum of M (typically 384 or 512) Fourier modes, 

M/2 

u{x,t)= J2 Cfc(i)e^^^ CfcGC. (45) 

k=-M/2+l 

Since u is real, we use the r2c version of the FFT algorithm, which only accesses the 
coefficients c^ with k > 0, assuming c_fc = c^. We also zero out the Nyquist frequency 0^/2 
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so that the total number of (real) degrees of freedom representing u at time t is M — 1. We 
use d = M/2 degrees of freedom to represent uq and T, namely 

Q = (ao,r,ai,&i,...,aM/4_i,6M/4-i) e I^'', {ck = ak + ibk, t = 0). (46) 

The remaining Fourier modes in uq are taken to be zero. The reason for using fewer Fourier 
modes in the initial condition is that in order to avoid ahasing errors, we want the upper 
half of the spectrum to remain close to zero throughout the calculation; therefore, we do 
not wish to give BFGS the opportunity to modify these coefficients. We increase M and 
repeat the calculation any time one of the high frequency (A; > M/4) Fourier modes of the 
optimal solution exceeds 10~^^ in magnitude at any timestep. 

To compute G{q), we write the Benjamin-Ono equation in the form 

ut = f{u) + g{u), g{u) = Hu,^, /(n) = - (^n^)^ , (47) 

where ■^u'^ is evaluated on the grid {xj = 2ttj/M : < j < M — 1} in physical space while 
H and dx are evaluated in Fourier space. The trapezoidal rule in physical space is used 
to evaluate the integral (j30p defining G. To evolve the solution, we use the stiffly stable, 
additive (i.e. implicit-explicit) Runge-Kutta method of Kennedy and Carpenter { KC031 
IWilOTb] known as ARK4(3)6L[2]SA with a fixed timestep h = T/N, where A^ is chosen to 
be large enough that further refinement does not improve the solution. Briefly, the idea of 
an ARK method is to treat / explicitly (as it is non-linear) while treating g implicitly (as 
it is the source of stiffness) : 



ki — J ( in + Cih, Un + h 2_^j O'ijkj + h 2_^j O-ij^j I 1 



^i = g \tn + Cih, Un + h Y,j a-ijkj + h ^j aij^.j j , 
Un+i = Un + h J2j bjkj + h J2j ijij- 



c 


A 


c 


A 




b^' 




6^ 


fo 


rf 


fo 


V g 



(48) 



The Butcher array for / satisfies aij = if i < j and for g satisfies aij = if i < j, 
which allows the stage derivatives to be solved for in order: ii, ki, £2, k2, ■ ■ ■ , ie, kQ, where 
our scheme has 6 stages. See [KC03] for the scheme coefficients and |Wil07bj for details 
on solving the implicit equations in the similar cases of Burgers' equation and the KdV 
equation. 

Once u{x, T) is known, we use the same scheme to solve the adjoint equation 

Ws = f{s,w)+g{w), g{w) = -Hwxx, f{s,w){x)=u{x,T-s)wx{x). (49) 

The main difficulty is that the intermediate stages of the ARK method require the value of 
u at intermediate times (between timesteps). For this we use cubic Hermite interpolation, 
matching u and ut at the timesteps straddling the required intermediate time: 

u{-,tn + eh) = (1 - 0)Un + eUn+1 " ^(1 " 9) [(1 - 29){Un+l " Un) " (1 " e)hdtUn + OhdtUn+l] 

where < < 1. This yields fourth order accurate values of u in the right hand side of 
l9|) . which is sufficient to achieve a fourth order accurate global solution w. We include 
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the option in our code to store u only at certain milemarker times, and then regenerate the 
data at ah timesteps between milemarkers as soon as the w equation enters that region; 
this dramaticahy reduces the memory requirements of the code at the expense of having to 
compute u twice. 

Once u{x,T) and w{x,T) are known with the period and initial conditions specified 
in g G M , we compute G{q) using the trapezoidal rule in physical space to evaluate the 
integral in (|30p . and we compute ^ by taking the FFT of -^ and scaling each component 



—-(x)ldx = 27r(^\ 



appropriately: 


OG 

dqo 


^271 

Jo 


dG 
dqi 


dG 
dT 


dG 


dG 


dq2k 


dak 



2tt 

[u{x, T) — no(x)]u((x, r) dx, 



27r 



use trap, rule 

in physical space 

(50) 



dG dG /"^^ 5G 



dq- 



2k+l 



w, - r i^(^) (-•'" - --") "- - ^' -{(£):} • <^ s ')■ 



We remark that these formulas for the derivatives of the numerical version of G essentially 
assume that we have solved the PDE exactly (so that the calculus of variations applies to 
our numerical solutions). This is reasonable in our case as we are using spectrally accurate 
schemes, but would cause difficulties if the numerical solution were only first or second order 
accurate in space or time. 

4.3 Choice of Penalty Function ip 

We still need to define the penalty function (p{uQ, T) in (j29p and show how to compute its 
gradient with respect to q. The purpose of ip is to pin down the mean and the phase shifts 
in space and time as well as to specify the bifurcation parameter. We have explored several 
successful variants which became more specialized as our understanding of the problem 
increased. As some of these variants may prove useful in other problems, we describe them 
here. 

Initially we did not include a penalty function in Gtot, but without it, the BFGS algo- 
rithm invariably converges to a constant solution. Next we constrained q2, the real part of 
the first Fourier mode ni(t) = ai{t) + ibi(t) at t = 0, to have a given value p. We reasoned 
that as long as p is not too large, the BFGS algorithm can vary q^ = bi (0) to find a periodic 
solution, so all we are doing is pinning down a phase. This was done by defining 

(^(no,r) = -([ao(0)-ao]^ + [ai(0)-/>]^) or ip{q) = -{[qo - ao? + [q2 - pf 

which works well to rule out the constant solutions but generally leads to traveling waves. By 
studying these traveling waves, we determined the formulas of Section [2] and also observed 
that for some choices of p and starting guess q^'^\ the wave becomes "wobbly," indicating 
that a non-trivial solution might be nearby. 



18 David M. Ambrose and Jon Wilkening 

To rule out traveling waves, we chose a parameter rj € [—1,1] and defined 
^{uo,T) = i([ao(0) -ao]2 + [ai(0) -p]2 + [ai(r/2) -7?ai(0)]2 + [6i(r/2) -7?6i(0)]2). 

Our idea here was that a (one-hump) travehng wave would have ry = ±1, depending on how 
many times it passed through the domain in time T; hence, intermediate values of r] would 
have to correspond to non-trivial solutions. To compute the gradient of if when it involves 
Fourier modes at later times, we simply solve another adjoint problem. Specifically, if if 
involves one of 

ak{T/2) = — / u{x,T/2)cos{kx)dx, bk{T/2) = — u{x,T/2)[-sm{kx)]dx, 

2vr Jo 27r Jq 

we will need to compute j^ak{T/2) or j^bf^{T/2), which can be done by setting 

wo{x) = —— cos{kx) , or wo{x) = — —— sm{kx) 
zvr zvr 

and solving ()4ip from s = to s = T/2; the result w{x,T/2) is the desired variational 

derivative. These may then be used to compute -^ak{T/2) or -^bk{T /2) as was done for 

G in (j50p . at which point it is a simple matter to obtain -^. 

This procedure proved very effective in obtaining non-trivial time periodic solutions. 
The BFGS algorithm is able to minimize Gtot down to 10~^^, at which point roundoff 
error prevents further reduction. With random initial data q^^' (which we tried before 
we had solved the eigenvalue problem for the linearization), the algorithm explores quite 
a wide region of the parameter space, with all components of q (including T) changing 
substantially — we do not seem to get stuck in non-zero local minima of Gtot- Once we do 
find a nontrivial solution, varying r] leads to other nearby periodic solutions. 

Studying this family of solutions, we finally realized that we were dealing with a four 
parameter family of nontrivial solutions with the mean, two phases and a bifurcation pa- 
rameter describing them. The main drawback of using -q as the bifurcation parameter is 
that the spatial and temporal phases are not specified independently, but instead depend 
on r/ in a complicated way. A more natural choice is to define 

^(no,r) = ^([ao(0) - ao]' + [afc(O) - pf + [6^(0)]^ + [dtak{Q)f), (51) 

i.e. we use if to impose the mean ao, the bifurcation parameter p, the spatial phase 6^(0) = 0, 
and the temporal phase dtak{0) = 0. Given any solution, we can always translate space and 
time to achieve the latter two conditions — we have not made any symmetry assumptions 
here. The index k we use depends on the number of humps N and bifurcation index n of 
the linearized solution; the only requirement is that zpf^k ^^ (|24p be non-zero. One readily 
checks that 

-I /'27r -| /'27r 

dtak{0) = — / Ut{x,0)coskxdx = — / \—k'^uo + {k/2)uQ]{—smkx)dx, (52) 
27r Jq 2tt Jq 

from which we obtain ^[9jafc(0)](x) = 2^(fc^ — kuo{x)) sinkx. Although ([^ does not rule 
out traveling waves, we have no difficulty bifurcating from traveling waves to non-trivial 
solutions by choosing a starting guess that includes first order corrections from the linear 
theory of Section [3l 
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5 Non-Trivial Time-Periodic Solutions 

We now use the methods described above to study the global behavior of non-trivial time- 
periodic solutions far beyond the realm of validity of the linearization about stationary and 
traveling waves. We find that these non-trivial solutions act as rungs in a ladder, connecting 
stationary and traveling solutions with different speeds and wavelengths by creating or 
annihilating oscillatory humps that grow or shrink in amplitude until they become part of 
the stationary or traveling wave on the other side of the rung. The dynamics of these non- 
trivial solutions are often very interesting, sometimes resembling a low amplitude traveling 
wave superimposed on a lower frequency carrier signal, and other times behaving like two 
bouncing solitons that repel each other to avoid coalescing. In this section, we present 
a detailed numerical study of the path of non-trivial solutions connecting the one-hump 
stationary solution to the two- hump traveling wave. In Section [6l we derive exact formulas 



for the solutions on this path. In a follow-up paper [AW], we classify all bifurcations from 
traveling waves, study the paths of non-trivial solutions connecting several of them, and 
propose a conjecture explaining how they all fit together. 

Consider the periodic solutions obtained by bifurcation from the one-hump stationary 
solution at the lowest frequency, wi^i. We arbitrarily set the mean oq = 0.544375 for these 
simulations (see Figure [l] above) , but as shown in Section [U any choice of ao < 1 would lead 
to similar results. In the top pane of Figure [3l we show the one- hump stationary solution 
lii^^(x) with /? = —\/{l — ao)/(3 — oq) together with the (initial conditions of the) two 
periodic solutions 

v'^°\x,t) = Re{zi,i(x)e^'^i*}, v'^^\x,t) = Im{zi^i{x)e''^''} (53) 



of the linearized equation (|T3j) corresponding to the first eigenvalue wi^i = 3 — oq of BA = 
—idx{Hdx — u). The natural period of these solutions is T = 27r/a;i^i = 2.55869. Note how 
the non-linearity of Benjamin-Ono distorts these two-hump perturbations as they travel (to 
the left) on top of the one-hump stationary "carrier" solution. Also note that v^^' and v^"^' 
are actually the same solution with a T/4 phase lag in time: 

t;(°)(x,r/4) = -v^^\x,Q) while v'^^^x.T/A) = u(°)(a;,0). 

We choose the real part of the first Fourier mode as the bifurcation parameter p so that 
/c = 1 in the definition ([^T]) of (p. As we vary p = ai(0) from — 2y^(l — qo)/(3 — ao) to 0, 
we traverse the trajectory from B to F in the bifurcation diagram of Figure HI The curves 
corresponding to the intermediate points HO, JO and KO along this path are shown in black 
in panes 2-4 of Figure [3l Along this path, we see that a second hump forms at x = 
while the center hump sharpens to accommodate the shorter wavelength of the two-hump 
traveling wave. If we instead increase \p\ near point B in the diagram, we obtain the lower 
path connecting B to G. Along this path, the center hump decreases in magnitude (curve 
H5), forms a dimple in the middle (curve JIO), splits into two humps (curve KIO), and 
again turns into a two- hump traveling wave (curve G). These curves are related to those 
on the path from B to F by a T/2 phase shift in time. If we change the sign of j3 (i.e. shift 
the phase by vr) in the stationary solution and call the resulting curve C, the bifurcation 
diagram is reflected about the T axis. The path from B (or C) to F is easier to compute 
due to the turning point in \p\ on the path from B (or C) to G. 
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Periodic solutions sampled at equal time intervals 




basis for periodic solutions 
of the linearized problem ■ 
T=2. 55869 



Figure 3: Progression from a one-hump stationary solution (top) to a two-hump travehng 
wave (bottom, moving left) by varying the real part of the first Fourier mode at t = while 
holding the mean oq constant and choosing the spatial and temporal phases such (/? in (j5ip 
is zero. The labels B, F, G, HO, JO, etc. correspond to Figures HHSl 
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Extrema of first Fourier mode vs. Period 
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Figure 4: Bifurcation from one-hump stationary solutions (B and C) to non-trivial time- 
periodic solutions that re-connect with two-hump traveling waves at F and G. 



Trajectories of first Fourier mode 




Figure 5: The trajectories of the first Fourier mode in the complex plane are exactly 
circular. The markers on the left lobe correspond to the solutions shown in Figure [3j 
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By the time we reach KO on the path from B to F, we can view our solution as a 
two-hump travehng wave with a small one-hump stationary perturbation corresponding to 
the first eigenvalue uj2,i = 1 in the linearization about the two-hump traveling wave. A full 
analysis of the linearization about traveling waves is given in the follow-up paper |AWj , but 
the idea is that if u{x) is a stationary solution and U{x,t) = u{x — ct) + c is a traveling wave, 
then the solutions v and V of the linearizations about u and U satisfy V{x, t) = v{x — ct, t). 
Now, the linearized solutions Re{z2,i(a^)e*'^^'^*} and Im{2:2,i(a;)e*'^^'^*} about the two-hump 
stationary solution have the property that Z2,i{x — vr) = —Z2,i{x)] hence, when they are 
used as perturbations on a two-hump traveling wave, they need to progress through an 
extra half-cycle in time to make up for the sign change. As a result, ij02^iT must belong to 
TT -|- 27rZ (rather than 27rZ itself) for the linearized solution to be periodic. It turns out that 
as we traverse the path from B to F, the period of the solution increases from T = 27r/a;i^i 
up to T = 7r/a;2,i = vr (rather than e.g. Svr or Svr). Note that as uj2^i = 1 is independent of 
ao, the path connecting the one hump stationary solution to the two- hump traveling wave 
always terminates at T = vr, regardless of the mean. 

In Figure O we plot the trajectories of the first Fourier mode ci(t) = ai{t) + ihi{t) in the 
complex plane for various choices of the bifurcation parameter p = ai(0). We were surprised 
to find that these trajectories are exactly circular; this will be discussed further below. The 
markers on the left (west) lobe of circles correspond to solutions plotted in Figure [31 for 
example, J19 corresponds to u (x, ^T), which is the dotted curve immediately to the right 
of the initial condition JO in the center pane of Figure [3l For visibility, we only plotted 10 
timeslices in the evolution of HO. 

The four parameter family of non-trivial solutions can be seen in Figure [5j A given 
solution is represented by one of the circular trajectories. The two main parameters de- 
scribing this family are the mean uq and the distance from the nearest point on the circle 
to the origin. A spatial phase shift of the initial condition by (with the sign convention of 
Eq. (j22p ) amounts to a clockwise rotation of the circle about the origin by 6 (or kO for the 
A;th Fourier mode). The north, east and south lobes of circles represent spatial phase shifts 
of the west lobe of solutions hy 9 = 7r/2, vr and — 7r/2, respectively, but any other phase 
shift E R is also allowed. Finally, a temporal phase shift amounts to choosing which point 
on the circle we assign to t = 0. Requiring that the initial condition have even symmetry 
yields either the west or east lobe of solutions with t = occurring along the real axis. 

We can also use other Fourier modes for the bifurcation parameter. This is especially 
important to track higher order bifurcations from multi-hump traveling waves — in these 
cases, the first several Fourier modes remain zero for all solutions in the family at all 
times t. But even for the simplest path connecting one- hump stationary solutions to two- 
hump traveling waves, it is useful to study other bifurcation diagrams representing this 
same family of solutions. In Figure El we show the result when the second Fourier mode 
is used instead of the first. By setting p = 02(0), we can now also see the bifurcation 
(labeled A) from the constant solution u = qq to the two-hump traveling waves; moreover, 
the points F and G that fell on top of each other in Figure Sj become distinct. The outer 
curve connecting F to G via A represents the set of two-hump traveling waves moving left 
with mean oq- This curve was plotted parametrically, setting 02 = ib2A^y^(l — a)/(3--a) 
and T = 2TT/[N{Na — oq)] with N = 2 and a ranging over all values such that a < 1 and 
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r<3.5. 

It is interesting to note that the bifurcation at F (and at G) from the two-hump travel- 
ing wave does not look like a pitchfork. Instead, the bifurcation curve enters at an oblique 
angle from one side only. This is because the second Fourier mode of the linearized solution 
v^^'{x,t) = Re{z2,i(x-|-t)e**} is zero (cf. (f2ll) above), so the first order effect on the bifurca- 
tion parameter p = 02(0) is zero as we move away from the two- hump traveling wave in the 
direction of v^^' . The derivative of T in this direction is also zero, so a heuristic argument 
based on the Liapunov-Schmidt reduction [GS85, IKie04j leads to an equation g{p, T) = 
for the bifurcation curve, where gf; = and ^^ = at the bifurcation. By contrast, the 
first Fourier coefficient of ^^'^•'(•,0) is non-zero and we do obtain a pitchfork bifurcation at 
F (and G) when we plot ai(0) vs. T, as was seen in Figure HI 

It turns out that the path of 02(0) vs. T from F to B is identical to the one from F to C; 
which one-hump stationary solution we end up with depends on whether we perturb the 
traveling wave in the direction of -\-v^^' or —v^^'. However, there is another direction we can 
move while keeping Gtot zero (with A; = 2 in (I5ip ). namely v^^'{x,t) = Im{2;2,i(x + t)e**}. 
This direction breaks the even symmetry of the initial condition, but the even Fourier 
modes still satisfy 6^(0) = and dtak{0) = 0; hence, the penalty function if does not 
exclude this direction when k = 2 in (I5ip . Depending on whether we perturb in the +f'^^ 
or —v^^' direction, we end up at either the one-hump stationary solution E, with maximum 
at X = 37r/2, or D, with maximum at re = 7r/2. This shows that our choice of penalty 
function (p in (j5ip does not rule out non-trivial solutions with asymmetric initial conditions: 
the solutions on the path from (F or G) to (D or E) are all asymmetric at t = 0; however, 
these solutions are related to the ones on the path from (F or G) to (B or C) by a phase 
shift in space. We have not found any periodic solutions that cannot be made symmetric 
at t = by such a phase shift. 

In Figure [71 we show the trajectories of the second Fourier mode in the complex plane. 
The markers labeled HO, JO, etc. again correspond to the solutions plotted in Figure [3j Un- 
like the first Fourier mode, these trajectories are not exactly circular — but by curve fitting 
we determined they are epitrochoids, (resembling Ptolemy's model of planetary motion, or 
"spirograph" trajectories), of the form 

C2(t) = C20 + C2,-ie'^' + C2,-2e*'"*, {00 = 27T/T) , (54) 

where the coefficients C2j (and oj) depend on the bifurcation parameter p. More generally, 
by curve fitting our numerical solutions, we have discovered a rather amazing property of 
solutions on this path: the /cth Fourier mode is found to be of the form 



Ck{t) = Y, Cfc.e-^J-*, {k>0, uj = 2tt/T), (55) 



where Ckj S M and C-k{t) = Ck{t). The general form of solutions on other paths connecting 
higher order bifurcations is similar, and is described in the follow-up paper [AWj. The four 
parameter family of non-trivial solutions is also nicely represented in this figure, where the 
parameters are the mean, the furthest point on the epitrochoid, a global rotation about the 
origin, and the choice of which point on the epitrochoid corresponds to t = 0. Note that 
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Figure 6: Bifurcation from the constant solution to a two-hump travehng wave and the 
path of non-trivial solutions connecting these to various one- hump stationary solutions. 
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Figure 7: The trajectories of the second Fourier mode in the complex plane are epitrochoids; 
see Equation ([5^ . The markers correspond to the solutions plotted in Figure [3l 
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a spatial phase shift of the initial condition by 6 leads to a rotation of a trajectory in this 
figure clockwise by 29, so the north and south lobes of circles in Figure [5] collapse onto the 
west family of epitrochoids (around D and E) in Figure [7] while the west and east lobes of 
Figure [5] collapse onto the east family here. 

6 Exact Solutions 

The discovery that the Fourier modes execute Ptolemaic orbits of the form (j55p led us to 
expect that it might be possible to write down the solution in closed form. In this section, 
we show how to do this for the path of non-trivial solutions connecting the one-hump 
stationary solution to the two-hump traveling wave. We have now learned of several other 
methods for finding exact solutions of Benjamin-Ono, notably the bilinear formalism used 
by Satsuma and Ishimori |SI79j and Matsuno [Ma t04j to construct multi-periodic solutions; 
the reduction by Case |Cas79| of the ODE ([9]) to a system shown by Moser to be completely 
integrable; and the approach of Dobrokhotov and Krichever |DK91| using the theory of finite 
zone integration to construct multi-phase solutions. Our approach highlights a feature of 
these solutions that has not been discussed previously, namely that the Fourier modes of 
these solutions turn out to be power sums of particle trajectories /3i{t) in the unit disk 
A C C whose elementary symmetric functions execute simple circular orbits in the complex 
plane. 

We start with the observation that the meromorphic solutions 

N 

u{x,t) = ao + ^n^,(t)(x), A(t) E A satisfying ([U]), (56) 

1=1 

have the property that the first A^ + 1 Fourier modes Ck{t) of u{x,t) are closely related to 
the trajectories of the /3;. Specifically, qq = cq is needed to write down the ODE ([TT]) . and 
we have 

/3i(t)+ • • • + /37v(t) = si{t), 2si{t) = ci(t), 

l3l{t)+---+l3J,{t) = S2{t), 2s2{t) = C2{t), 

(57) 

/3f (0+ • • • + t^Nit) = SN{t), 2sN{t) = CN{t). 

It is a standard theorem of algebra [ vdW70] that the elementary symmetric functions 

^j= E Ai---A,. (j = l,...,iV) (58) 

h<---<lj 

are polynomials in the power sums, e.g. 

, sf - S2 si- 3siS2 + 2S3 
0-0 = 1, (Tl=Si, CT2 = , (T3 = . (59) 



The general recurrence relation is 



o-Q 



-^{-ly-'Tj-iSi, iJ = l,...,N). (60) 



^1=1 
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Figure 8: Left: trajectories of /3i, /32 for the solutions in Figures [3HZ1 As we vary the 
bifurcation parameter, the trajectories change from two disjoint, counterclockwise loops to 
one larger orbit in which /?i, /32 exchange positions over the course of one period. Right: 
the trajectories of ci, (T2 are exactly circular. 



The /3/ are then the zeros of the polynomial 

AT N 



(61) 



1=1 



j=0 



We can test whether a given numerical solution u{x,t) is an A^-particle meromorphic solu- 
tion by computing its first A'' -|- 1 Fourier coefficients Cfc(O) = 2sfc(0), using (i60]l to obtain 
the symmetric functions o"j(0), solving for the roots A(0) of the polynomial on the right 
hand side of (j6ip . and checking that higher power sums do in fact agree with the Fourier 
coefficients of the solution: 



/3i(0) + --- + /3^(0) = -Cfc(0), {k>N + l) 



(62) 



Using this approach, we find (numerically) that the solutions on the path connecting 
the one- hump stationary solution to the two- hump traveling wave are 2-particle solutions. 
Moreover, the trajectories of the first two symmetric functions appear to be of the form 

ai=pi+l32 = -A + Be''^\ (63) 

CT2 = A/32 = -Ce^"*, (64) 

where A, B, C, oo are positive constants; see Figure [HI We now prove this rigorously. 

Theorem 1 There is a four-parameter family of time-periodic, two-particle solutions of the 
form, 

u{x,t) = ao + Ui3^^t){x) + up^^t){x), (65) 
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where /3i(t) and /?2(i) are the roots of the equation 

z^-ai{t)z + a2it) = (66) 



and 



cji(t) = [-A + 5e^'^(*-*«)]e-^^ a2it) = [-Ce^'^(*-*o)]e-2ie^ (g^^ 

(3 - ao)^[(3 - ao) - (7 - ao)C^] [(1 - ao) - (5 - ao)C^] 
^ ^ (3-ao)2-(5-ao)2C2 ' ^^^^ 

i? = I^^AC, (69) 

6 — ttQ 



(3 - ao)' - (5 - aoYC 



(70) 



(3 - ao) - (5 - Qo)C2 • 

T/ie four parameters are the mean ao < 1, iw^o phases 6,tQ € M, and a real number C 
ranging from C = (at the one-hump stationary solution) to C = J ^~^° fai the two-hump 
traveling wave). 

Proof: It suffices to consider the case tliat ^ = and to = as the general case fohows 
immediately. If we try to substitute /3i^2 = ^ ± 2 \/^i ~ ^^"^ '™-'^^ ^^^ system 

/3i = ^o ^ + -^ ^=i+-R 5^ + <3-ao)/3i (71) 

Pi - P2 Pi - Pi /?i - P2 

■ -2i/3| 2i/3| 2i/?f , ,n 

/32 = ^5-^ + -^^ + ^^ + <3-ao)/32 (72) 

P2 - Pi P2 - Pi P2 - P2 

and solve for A, B and uj in terms of C and ao, the algebra becomes unmanageable. However, 
we can re-write this system in terms of ai and o"2 to obtain 

..[ /9i(/3iA) , /3i(/3i^2) , /32(/32^i) , /32(/32^2) l , .,, . 

(Ti = — 2i < = = = ^- > + i 1 — an Cn , 

ll-/3i/3i l-/3i/32 l-/32/3i 1-/32/32 J ^ ^ ' 

^.f /3iA , /3i/g2 , P2h , /32^2 1 ^^... . 

0"2 = — 2z < =- -\ = = ^- > 1T2 + 2i 2 — an]CJ2. 

1 1 - /3i/3i ^ 1 - /3i/32 ^ 1 - /32/3i ^ 1 - /32/32 / ' ^ ^ ^^ ^ 

The expressions inside braces remain invariant if we interchange /3i and /32; hence, they 
may be written as rational functions of ui, (J2, (^i, ^2- Explicitly, we have 

-Pi -P2 

(Ti = -2i— +i(l -ao)cri, 0-2 = -2i— 0-2 + 2i(2 - ao)cr2, (73) 

Pi = CTiO-i - 2a-icr2 - 2a\a2 + 6cri|(T2p - (Jia\a2 + 2a\ai\a2\^ - 2aia2d-2 - 2ai\a2\'^, 
P2 = ki|2(l + 3|cJ2p) + 4|cj2p(l - |a2p) - 2{al(T2 + aja2), 
Q={1- \a2ff - |ai|2(l + \a2f) + {ala2 + ^1^2). 
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Since Q is a product of non-zero terms of the form (1 — Pil3j), it is never zero. If we assume 
ai = -A + Be''^^ a2 = -Ce*^*, and C / 0, we find that ^ holds as long as 



2Pi + (1 - ao)^iQ +^ - 2P2^2 + (4 - 2ao)a2Q =0, (74) 



CI 



-2P2 + (4 - 2ao - uj)Q = 0. (75) 



We eliminated lo in ()74p using ai = icoBe^^^ = —^ct2- Next, we collect terms containing 
like powers of e*'^* and set them each to zero. This yields 7 polynomial equations in the 
variables A, B, C, ao and m; however, several of them are redundant due to relationships 
such as Q(-^) = Q(^) in the decomposition Q = Q^-^'^e''^^ + Q(°) + Q^^^e*'^*. Equation 
(|74p yields 4 such equations; two of them are satisfied if we choose B as in ()69p while the 
remaining two are satisfied if we also choose A as in (j68p . With these choices, all three 
equations associated with ([75|) are satisfied provided w satisfies ([70|) . The special cases 

{C = 0, A = Jj^, B = 0} and {C = J^E^, A = 0, B = 0} are seen to correspond to 
the one-hump stationary solution and two-hump traveling wave, respectively, as discussed 
in Section [2j D 

We have verified that the curve connecting B to F in the bifurcation diagram of Figure S] 
is recovered if we set oq = .544375 and plot 2{—A + B) versus T = 2tt/u} using the above 

formulas for A, B and uj with C ranging from to a / g~^° . 



7 Conclusion 

We have presented a general method for finding continua of time-periodic solutions for 
nonlinear systems of partial differential equations. We have used our method to study global 
paths of non-trivial time-periodic solutions connecting stationary and traveling waves of the 
Benjamin-Ono equation. In spite of the non-linearity and non-locality of the Benjamin-Ono 
equation, these non-trivial solutions can be interpreted as distorted superpositions of the 
stationary or traveling waves at each end of the path. Our numerical method is accurate 
enough that we are able to use data fitting techniques to recognize the analytical form of 
the solutions. In particular, the Fourier coefficients Ck{t) of these solutions follow Ptolemaic 
orbits of the form (j55p . This led us to reformulate the equations governing meromorphic 
pole dynamics to reveal an exact formula for the solutions on the four-parameter path 
connecting the one-hump stationary solution to the two-hump traveling wave. 

In the future, we plan to apply this method to more complicated systems arising in fiuid 
dynamics, namely the vortex sheet and water wave problems. This will allow for comparison 
with prior numerical and analytical results |HLS97) . [PTOlj . |IPT05] . Additionally, as the 
Benjamin-Ono equation is meant as a model for internal waves in a deep, stratified fiuid, 
it will be of interest to compare time-periodic vortex sheets and water waves with time- 
periodic solutions of Benjamin-Ono. 
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